# Esquema de Diferenças Finitas para Fluxo Laminar;
# Raios nas extremidades da cela computacional;
# Velocidades no centro da cela computacional;


dpdz = -1;
R = 1;

m = 100;

dr = R/(m+1);
dw = dr;

for i = 1:(m+1)

	mu(i) = 1;
	r(i) = (i-1)*dr;
	h(i) = dw;

end

# Lado Direito da Discretização

LD = zeros(m,1);

for i = 1:m

	LD(i,1) = (dpdz)*r(i);

end

# Lado Esquerdo da Discretização

LE = zeros(m,m);

dw1 = r(2) - r(1);

coef_11 = -mu(1)/(dw1);
coef_1mais = mu(2)/(dw1);

dwm = r(m)-r(m-1);

coef_mi = -2*r(m)*mu(m)/(dwm*dwm);
coef_mmenos = -(mu(m) + (mu(m)-mu(m-1))*r(m)/(2*dr))/(2*dwm) + mu(m)*r(m)/(dwm*dwm);

LE(1,1) = coef_11;
LE(1,2) = coef_1mais;
LE(m,m) = coef_mi;
LE(m,m-1) = coef_mmenos; 

for i = 2:(m - 1)

	coef_ii = -2*r(i)*mu(i)/(h(i)*h(i));
	coef_imenos = -(mu(i) + (mu(i)-mu(i-1))*r(i)/(2*dr))/(2*h(i)) + mu(i)*r(i)/(h(i)*h(i));
	coef_imais = (mu(i) + (mu(i)-mu(i-1))*r(i)/(2*dr))/(2*h(i)) + mu(i)*r(i)/(h(i)*h(i));

	for j = 1:m

		if j == i

		LE(i,j) = coef_ii;

		end
		
		if j == (i - 1)

		LE(i,j) = coef_imenos;

		end

		if j == (i + 1)

		LE(i,j) = coef_imais;

		end

	end

end

W = LE\LD;

max_v = 0;

for i = 1:m

	if W(i) >= max_v
		
		max_v = W(i);

	end

end

